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ABSTRACT 

Progress  has  been  made  in  both  algebrsuc  and  elliptic  grid  generation,  (i)  The  control 
point  form  (CPF)  of  algebraic  grid  generation  has  been  improved  in  two  aspects.  First, 
new  blending  functions  are  developed  which  allow  a  computer  software  user  to  specify  the 
locations  of  control  surfaces  arbitrarily.  Second,  a  new  implementation  scheme  is  used  in 
the  CPF  to  recapture  the  clustering  feature  of  given  grids,  (ii)  An  effective  technique  of 
curvature  control  in  elliptic  grid  generation  has  been  developed  and  implemented  to  control 
the  grid  point  distribution  along  curved  boundaries. 
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CHAPTER  I.  INTRODUCTION 


The  control  point  form  (CPF)  of  algebraic  grid  generation  for  a  2D  (or  3D)  grid 
results  from  a  combination  of  the  multisurface  transformation  [1,2]  and  the  transfinite 
interpolation  [3].  The  multisurface  transformation  can  be  thought  of  as  a  new  kind  of 
curve  generation  technique,  hke  many  other  kinds  of  curve  generation  methods.  From  the 
point  of  view  of  a  computer  software  user,  the  user  would  like  to  specify  the  locations 
of  control  points  at  his  or  her  choice.  This  fleability  for  a  user  requires  that  a  set  of 
blending  functions  can  be  constructed  for  any  given  set  of  control  points.  Previously,  the 
blending  functions  used  in  the  multisurface  transformations  have  been  developed  [1,2]. 
However,  the  locations  of  control  points  in  the  physical  space  are  determined  after  the 
blending  functions  are  constructed.  As  a  result,  the  locations  of  the  control  points  in 
the  physical  space  cannot  be  chosen  arbitrarily.  Put  the  problem  into  the  context  of  grid 
generation  within  a  2D  physical  region,  certainly  we  can  provide  a  set  of  control  points  for 
a  user.  From  the  point  of  view  of  a  computer  software  user,  he  or  she  may  still  want  to 
insert  or  remove  one  or  more  control  curves.  In  order  to  meet  the  user’s  needs,  we  should 
develop  a  new  set  of  blending  functions  which  would  allow  a  user  to  choose  the  locations 
of  control  curves  arbitrarily.  Once  we  have  such  a  set  of  blending  functions,  we  can  easily 
handle  the  problem  of  allowing  a  user  to  add  and  delete  one  or  more  control  curves.  In 
Chapter  II  of  this  report,  we  present  our  results  on  new  blending  functions-both  C^-  and 
C^-continuity-which  allow  arbitrary  specification  of  the  locations  of  the  control  points. 

The  blending  functions  we  mentioned  in  the  paragraph  above  are  all  local  interpolation 
functions.  This  locality  gives  the  CPF  the  capability  of  changing  or  modifying  a  given  grid 
locally.  The  general  procedure  goes  like  this:  A  user  has  a  grid  which  has  been  generated 
in  one  of  other  methods.  Then,  the  user  wants  to  improve  the  quality  of  his  grid  locally 
using  the  CPF  method.  Usually,  this  requires  first  regenerating  the  user’s  grid  using  the 
CPF  method.  After  regeneration,  the  user  can  move  one  or  more  control  points  to  change 
the  grid  locally.  In  previous  application  of  the  CPF  to  regenerating  a  given  grid,  it  was 
found  that,  for  a  given  grid  with  clustering  along  one  or  more  boundary  curves,  the  grid 
regenerated  using  the  CPF  method  often  fails  to  recapture  the  general  clustering  feature  of 
the  given  grid  near  concave  regions  (cf.  Figs.  3. 1-3.3).  To  make  the  application  of  the  CPF 
method  successful,  we  need  first  to  find  out  the  reason  of  the  failure  and  then  to  correct  it. 
It  turns  out  that  the  reason  of  the  failure  was  in  the  usual  implementation  of  connecting  the 
continuous  world  of  mathematics  (analytic  formula)  and  the  discrete  world  of  numerical 
computations  (usually  encountered  in  numerical  grid  generation).  For  this  problem,  we 
do  not  need  any  new  blending  functions.  For  any  set  of  blending  functions,  whether  old 
or  nev:,  we  can  use  them  to  recapture  the  gener^J  clustering  pattern  of  any  given  grid 
along  one  or  more  boundary  curves,  since  they  all  satisfy  the  uniformity  condition  [l].  In 
Chapter  III  we  report  our  progress  made  in  this  area  and  our  new  implementation  of  the 
control  point  form  of  algebraic  grid  generation.  Also  in  Chapter  III,  the  grids  generated 
using  both  old  and  new  implementation  methods  of  the  CPF  are  presented.  From  those 
illustrations  of  2D  grids,  one  can  easily  .soe  th**  dramatic  improvements  in  the  concave 
regions  for  a  grid  with  clustering. 

Chapter  IV  of  this  report  deals  with  improving  the  quality  of  grids  generated  ellipti- 
cally.  It  describes  an  effective  technique  for  controlling  the  tendency  of  elliptic  systems  to 
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pull  grid  points  away  from  concave  regions  of  the  grid.  This  tendency  is  responsible  for 
the  poor  grid  resolution  near  concave  regions,  exhibited  by  elliptic  grids.  The  discussion  in 
Chapter  IV  takes  a  close  and  detailed  look  at  the  discrete  form  of  the  governing  equations 
and  isolates  the  terms  which  pull  the  points  in  various  directions.  In  other  words,  Chapter 
IV  describes  a  technique  for  improving  grid  quahty  by  controlling  the  effects  of  curvature 
of  the  boundaries.  The  term  used  to  designate  this  technique  in  this  report  is  “curvature 
control.”  Results  from  a  program  which  implements  the  above  technique  are  presented 
and  show  improvements  in  grid  quality  over  grids  generated  without  the  help  of  curvature 
control. 
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CHAPTER  II.  MULTISURFACE  TRANSFORMATION 

Ill  this  chapter,  we  report  our  new  hlendiiig  functions  developed  in  the  Phase  II  work. 
This  chapter  is  organized  in  the  following  way.  In  Section  2.1,  we  give  the  general  properties 
of  the  multisurface  transformation.  In  Section  2.2  we  review  previous  results  on  C^- 
continuity  blending  function.  In  Section  2.3  we  construct  new  the  C^-continuity  blending 
functions.  In  Section  2.4  we  present  old  results  on  C^-continuity  blending  functions.  In 
Section  2.5  we  develop  new  C^-continuity  blending  functions. 

2.1.  General  Formalism 

To  state  the  results  mathematically  some  notation  is  needed.  For  this  purpose,  let  Pj, 
P2,  ....  Pat  be  the  given  sequence  of  points  in  a  2D  space;  let  r  be  the  curve  parametrization; 
let  P(r)  be  the  position  at  r  along  the  desired  curve;  let  rj,  r2,  ...  be  the  successive 

parametric  locations  to  interpolate  the  directions  of  (P2  —  Pi),  (P3  —  P2),  •••,  (Pn  — 
Pa?^_i);  and  let  V’i(^),  V’2(^),  •••,  be  the  corresponding  interpolation  functions 

which  successively  separate  each  direction  by  assuming  a  non-zero  value  at  the  associated 
location  while  vanishing  at  the  remaining  locations  for  interpolation.  In  two  dimensions 
Pfc  =  (xk-iVk)  and  P(r)  =  (i(r),t/(r)).  With  this  notation  the  curve  is  given  by 

where 

=  f  '4>ki^)dx,  k  =  1,2,  ...,N  —  1.  (2.1.2) 

Jri 

To  witness  the  basic  specifications  mentioned  above,  it  is  easy  to  check  the  end  conditions 
P(ri)  =  Pi  and  P(rAr-i)  =  Pyv  and  the  interpolatory  condition  that  dP{rk)/dr  is  in  the 
direction  of(Pfe+i— Pfc)  for  each  k  from  1  to  iV  — 1.  In  the  context  of  coordinate  generation 
for  two-  or  three-dimensional  regions,  the  endpoints  become  boundary  surfaces  and  the 
interior  points  becomes  control  surfaces.  For  this  reason  the  transformation  generated  by 
curves  of  the  above  form  has  been  called  a  multisurface  transformation. 

Corresponding  to  the  N  points  Pi,  P2,  ...,  Pat,  we  label  their  parametric  values  as 
61,  62,  ....  67V,  which  are  in  an  increasing  sequence 

61  <  62  <  ...  <  b^.  (2.1.3) 

The  values  of  hi,  62,  •••,  b^  can  be  specified  or  determined  by,  say,  bk  =  Pfc  •  t  (we  will 
discuss  how  to  choose  these  6’s  values  in  Chapter  III).  Because  of  relation  (2.1.3),  the  TV  —  1 
parametric  values  for  TV  —  1  derivative  directions  are  also  in  an  increasing  sequence, 

Ti  <r2  <  ...  <  Tf^-i.  (2.1.4) 

The  two  sets  of  parametric  values  are  coincide  at  two  ends  of  the  curve, 

bi=Ti,  bN  =  rN-i,  (2.1.5) 
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The  interpolation  function  i>k{r)  in  Eq.  (2.1.2)  satisfies  a  cardinality  condition, 


ipk{rj)  =  Sjk,  j,k  =  1,2,..., N  -  1.  (2.1.6) 

An  important  condition  for  the  curve  (2.1.1)  is  a  uniformity  condition,  which  states 
that,  when  projected  on  to  a  vector  t,  the  curve  becomes  a  simple  linear  curve  in  the 
parametric  space,  P(r)  -  t  =  r.  Let 


Ck  =  bk+i  —  bk,  k  =  1, 2, ...,  N  —  1, 
the  uniformity  condition  is  expressed  mathematically  as 

N-i 


=  Ti  +  Yl 


Gk{r) 


Ck. 


^  GkirN-i) 

Taking  derivative  with  respect  to  r  on  both  sides  of  Eq.  (2.1.8),  we  obtain 

N-l 


(2.1.7) 


(2.1.8) 


^k(r) 


Ck. 


(2.1.9) 


“  Gkirtf-i) 

The  cardinality  condition  (2.1.6)  makes  the  summation  collapse  to  only  one  term,  yielding 


Ck  =  Gk{rN-i),  ib  =  l,2,...,jV-l.  (2.1.10) 

Substituting  Eq.  (2.1.10)  back  into  Eqs.  (2.1.8)  and  (2.1.9),  we  simplify  the  uniformity 
condition  to 


N-l 


r  =  ri  +  ^  Gkir), 
k=l 

(2.1.11) 

N-l 

fc=i 

(2.1.12) 

The  actual  values  of  the  parameters  and  bk  are  related  in  some  way.  If  one  first  chooses 
Tfc’s,  then  bk's  are  determined  by  Eqs.  (2.1.7)  and  (2.1.10).  On  the  other  hand,  if  one  first 
is  given  bk's,  then  one  has  to  decide  rfc’s  by  using  Eq.  (2.1.10)  and  other  choices.  We  will 
discuss  these  “other  choices”  in  Sections  2.3  and  2.5.  In  any  case,  since  the  two  sets  of 
parametric  locations  are  coincide  in  two  ends  [see  Eq.  (2.1.5)],  the  sum  of  all  intervals  in 
terms  of  one  set  of  parametric  locations  must  be  equal  to  that  in  terms  of  the  other  set  of 
parametric  locations.  In  other  words,  if  we  define 

hk  =  ^fc+i  “  fki  —  l,2,...,iV  —  2,  (2.1.13) 

then  we  have  a  restriction  for  the  two  sets  of  intervals. 
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In  the  rest  of  this  chapter,  we  report  C®  and  C*  interpolants  V’jt,  i.e.,  and  C~ 
blending  functions  Gfc.  We  first  consider  continuity  for  blending  functions  Gfc(r)’s  in 
Sections  2.2  and  2.3:  The  interpolation  functions  i/>fc(r)’s  are  piecewise  continuous  them¬ 
selves,  whereas  the  blending  functions  Gfc(r)’s  arc  continuous  up  through  first  derivatives. 
Then,  we  consider  blending  functions  in  Sections  2.4  and  2.5.  The  interpolation  func¬ 
tions  ipki^Ys  are  continuous  up  through  first  derivatives,  whereas  the  blending  functions 
Cjfc(r)’s  are  continuous  up  through  second  derivatives. 

2.2.  Previous  Results  on  C^-Continuity  Blending  Functions 

The  simplest  local  interpolants  belong  to  the  class  C®  of  continuous  functions.  The 
interpolation  functions  ipk  consists  of  two  pieces  for  k  =  2,3,...,W  —  2.  Each  of  the 
interpolants  V"!  and  V’AT-i  for  the  two  end  points  is  made  up  of  only  one  piece.  An 
illustration  of  continuous  piecewise  linear  interpolation  functions  is  plotted  in  Fig.  2.1. 


Figure  2.1.  Continuous  piecewise  linear  interpolation  functions. 

In  the  region  r/t  <  r  <  r^+i,  only  two  blending  functions  are  non-zero,  and  the 
uniformity  condition  (2.1.12)  reduces  to 


+  V’fc+i(T’)  =  1,  rk<T<rk+i,  k  =  -2.  (2.2.1) 

For  arbitrary  positions  of  r*,  the  interpolation  functions  must  be  of  the  following  form, 


V’i(»*) 


/i(xi),  forr,  <r<r2, 
for  r2  <  r  < 


V’fc(T')  = 


0,  for  ri  <  r  <  rt-i, 

fk-i{xk-i),  for  Tk-i  <T  <rk, 
I  -  fkixk),  for  Tfc  <  r  <  rjk+i, 
0,  for  rfc+i  <  r  < 


for  k  =  2,3,...,  AT  -  2, 


for  ri  <  r  <  r/^_2, 
for  <r  <  rs~i, 


(2.2.2a) 


(2.2.2b) 


(2.2.2c) 
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where 


T  —  Tk 

Xk  =  - 

’"fc+l  —  T'k 


r  -  Tk 
hk  ' 


k  =  -  2. 


(2.2.3) 


In  Ref.  [l],  simple  piecewise  linear  interpolation  functions  were  used  (see  Fig.  2.1), 


fkix)  =  X.  (2.2.4) 

Substituting  Eqs.  (2.2.2)  and  (2.2.4)  into  Eq.  (2.1.2)  and  completing  the  integr2ds,  we 
obtain  the  blending  functions 


{hi{xi-\x\),  forri<r<r2, 
\\hi,  for  ra  <  r  <  rjv_i, 


Git(r) 


'0, 

2  ®  fc  — 1  ^ 

^hk-1  +  hk{xk  -  \xl), 
^  2{^k-l  +  hk). 


for  Ti  <  r  <  Tk-i , 
for  Tk-i  <r  <rk, 
for  Tk  <r  <  Tk+i, 
for  rjk+i  <r<  rs-i, 


for  k  —  2,3, ...,  N  —  2,  and 


(2.2.5a) 


(2.2.5b) 


GN-i{r) 


Jo,  for  ri  <  r  <  ts-2^ 

\\hN-2xlf^2^  for  rjv-2  <  r  <  rA/_i. 


(2.2.5c) 


A  graphical  illustration  of  the  blending  functions  Gfc(r)  is  drawn  in  Fig.  2.2.  Each  function 
increases  monotonically  from  0  to  its  maximum  value  Gk{TN-i ).  According  to  Eq.  (2.1.10), 
these  maximum  v^dues  give  the  intervals  of  the  control  points  Pfc  in  the  parametric  space, 


Cl  = 


(2.2.6a) 


Ck  =  -{hk-i  +  hk). 


k  =  2,3,...,N  -  2, 


(2.2.6b) 


Cj^-i  =  -h7v-2-  (2.2.6c) 

Since  the  coefficients  represented  by  Eqs.  (2.2.5)  are  known  functions,  the  curve  of  equation 
(2.1.1)  depends  only  upon  the  sequence  of  points  Pi,  Pa,  ...,  Pn- 

The  values  of  rt’s  can  be  chosen  arbitrarily.  The  simplest  choice  is  to  choose  a  unit 
spacing  in  the  interpolation  points,  Tk  =  k.  This  causes  both  the  interpolation  functions 
and  their  integrals  Gk  to  be  translations  of  a  function  about  the  origin.  The  overall 
consequence  is  a  simply  stated  curve  definition.  In  an  analytical  form  the  origin- centered 
interpolation  function  is  given  by 


V»(x 


1  - 
0, 


for  — 1  <  X  <  1, 
otherwise  , 


(2.2.7) 
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and  its  integral  is  given  by 


i7(i) 


0,  for  I  <  —  1, 

|(i  +  1)^,  for  — 1  <  X  <  0, 
1  —  \{x  —  1)^,  for  0  <  I  <  1, 

1,  for  X  >  1. 


In  this  simple  case,  the  coefficients  in  Eq.  (2.1.1)  are  given  by 


Gi{r) 

Gi{rf4-\) 


=  2fl(r  -  1)  -  1, 


Gkjr) 

Gfc(rjv-i) 


=  S7(r  —  fc). 


k  =  2,3,...,^■  -  2, 


Gn-\{t') 

Gn-i{tn-i) 


=  2il(r-iV  +  l). 


Also,  the  intervals  Ck  are  of  the  form  of  “half  one  one  ...  one  half,” 


Cl 


=  Gn~2  = 


1 

2’ 


(2.2.8) 


(2.2.9a) 


(2.2.9b) 


(2.2.9c) 


(2.2.1Ca) 


Cfc  =  l,  fc  =  2, 3, ...,  AT  -  2.  (2.2.10b) 

Besides  piecewise  linear  interpolation  function  (2.2.4)  or  (2.2.7),  it  is  certainly  true  that 
one  can  construct  various  other  C°  continuity  interpolation  functions.  For  example,  one 
can  have  the  following  piecewise  trigonometric  interpolation  function, 


for  — 1  <  a:  <  1, 
otherwise  , 


(2.2.11) 


which  is  continuous  up  through  first  derivatives.  Thus,  the  corresponding  blending  func¬ 
tions  Gfc(r)  are  continuous  up  through  second  derivatives. 

The  scheme  discussed  in  this  section  determines  the  two  sets  of  parametric  locations 
{rk,hk)  and  (6t,C*)  in  the  following  order; 


{‘>’k,hk)—^{ikk{'r),  Gjt(r))— >(6fc,  Cfc).  (2.2.12) 

From  the  viewpoint  of  application,  this  kind  of  scheme  has  a  disadvantage:  If  a  computer 
software  user  wants  to  specify  the  parametric  locations  bk  at  his  or  her  choice,  then  it  does 
not  tell  us  how  to  determine  rjt  and  Gfc.  Besides,  it  is  not  possible  to  find  a  corresponding 
set  of  Tfc’s  in  certain  situations.  For  example,  since  hk  >  0  for  all  k,  we  see  from  Eqs. 
(2.2.6)  that  Gi  <  C2  always.  Thus,  this  scheme  is  not  capable  of  handling  the  case  in 
which  Cl  >  G2.  More  flexible  scheme  is  needed.  We  will  present  a  new  sc’  ome  in  the  next 
section. 
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2.3.  New  Results  on  C^-Continuity  Blending  Functions 

For  any  given  set  of  fcjt’s,  the  question  is  how  to  determine  rjt's.  Our  way  of  choosing 
rjt’s  is 


hi  =  Cl  +  -(72, 

(2.3.1a) 

hk  = -(Ck  +  Ck+i),  fc  =  2,3,  ...,iV  -  3, 

(2.3.1b) 

hN-2  = -Cn-2  +  Cn-1- 

To  consol'date  Eqs.  (2.3.1)  into  a  single  equation,  we  introduce 

(2.3.1c) 

Ck  =  -Ck,  fc  =  2,3,  ...,iV  —  2, 

(2.3.2a) 

Ck  =  Ck,  k  =  l,N-l, 

and  thus  have  a  simple  form  for  Eqs.  (2.3.1), 

(2.3.2b) 

hk  =  Ck  +  Ck+i,  fc  =  1,2,  ...,iV  -  2. 

(2.3.3) 

The  piecewise  linear  function  described  in  Section  2.2  does  not  work  here,  since 
Gk{rN-i)  =  ^{hk-i  +  hk)  in  Eq.  (2.2.3b)  does  not  satisfy  Eq.  (2.1.10)  in  general.  More 
general  and  flexible  interpolation  functions  are  needed.  Starting  from  the  first  interval 
rj  <  r  <  r2  for  the  interpolation  functions,  according  to  Eqs.  (2.1.2)  and  (2.1.10),  we  have 
to  require  that  the  area  under  the  interpolation  function  is 


i>i(r)dT  =  Gi(r2)  =  Gi{rN-i)  =  Ci  =  Cy. 


(2.3.4) 


Because  of  Eqs.  (2.3.4)  and  (2.3.1a),  this  makes  the  area  under  the  interpolation  function 
V’2(^)  within  the  first  interval  is 


^2(^2) 


')dr= 

J  ri 


^i{r)]dr  =  hi  —  Cl  =  -C2  =  C2. 


(2.3.5) 


Then,  we  examine  the  second  interval  r2  <  r  <  rz.  In  order  to  satisfy  Eq.  (2.1.10),  we 
have  to  require  that  the  area  under  the  interpolation  function  V’2(^)  within  the  second 
interval  is 


J  ^2{r)dr  =  Gzirz)  —  G2(»’2)  =  ^zir^-i)  —  -(72  =  -C2  =  C2. 


(2.3.6) 
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Proceeding  in  this  way,  we  find  that,  for  the  kih  (k  —  l,2,...,iV  —  2)  interval,  the  require¬ 
ments  are 


Gk{rk+i)  —  <5fc(rjt)  = 


ifk{r)dr  =  Cfc, 


(2.3.7a) 


Gk+ii^k+i)  =  I  '^k+iir)dr  =  Ck+i-  (2.3.7b) 

Jr* 

Equation  (2.3.7b)  can  be  rewritten  as 

Gkirk)  =  Ck,  ik  =  2,3,...,JV  - 1.  (2.3.8) 

However,  Gi(ri)  =  0  is  an  exception. 

We  have  mentioned  above  that  the  piecewise  linear  interpolation  function  will  not 
satisfy  Eqs.  (2.3.7)  in  general.  In  order  to  satisfy  Eqs.  (2.3.7),  one  more  parameter 
(freedom)  should  be  allowed  in  the  form  of  the  interpolation  function.  One  choice  is  to  use 
a  quadratic  interpolation  function, 


fk{x)  =  X  +  Akx{x  -  1).  (2.3.9) 

Equation  (2.3.9)  has  been  set  in  a  convenient  form  which  satisfies  the  two  requirements 
ipkirk)  =  1  —  /fe(0)  =  1  and  tkkirk+i)  =  1  —  /fc(l)  =  0-  Substituting  Eqs.  (2.2.2)  and 
(2.3.9)  into  Eqs.  (2.3.7),  it  is  straight  forward  to  find  that 

Ak  =  3iCk  -  Ck+i)/hk,  =  1, 2, ...,  N -2.  (2.3.10) 

The  sign  of  Ak  depends  on  whether  Ck  >  Ck+i  or  Ck  <  Ck+i-  When  Ck  =  Ck+i,  the 
coefficient  Ak  vanishes,  Ak  =  0.  The  quadratic  interpolation  function  fk{x)  is  not  positive 
definite  within  the  range  0  <  x  <  1.  For  example,  when  Ak  >  1,  i.e.,  when  Ck  >  2Ck+i, 
the  function  fk  =  —{Ak  —  l)^ /4Ak  is  negative  at  x  =  {Ak  —  l)f2Ak  <  To  ensure  the 
interpolation  functions  to  be  positive  definite,  we  can  choose  another  form 


/fc(x)  =  x”**,  (2.3.11) 

which  has  the  same  end  values  /fc(0)  =  0  and  /*(!)  =  1  as  the  /fc(x)  given  in  Eq.  (2.3.9). 
Substituting  Eqs.  (2.2.2)  and  (2.3.11)  into  Eqs.  (2.3.7),  we  obtain 

mk  =  -^,  fc  =  l,2,...,W-2.  (2.3.12) 

Gk+i 


Integrating  the  interpolation  functions  ^fc(^)?  we  obtain  the  blending  functions 

=  ^  ^  ^ (2.3.13a) 

(Cl,  for  rj  <  r  <  tn-i, 
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Gfc(r)  = 


(2.3.13b) 


0,  for  ri  <  r  <  rt-i, 

Hk-i{xk-i),  for  Tk-i  <r  <rk, 

Ck  +  hkXk  -  Hkixk),  for  Tk  <r  <  Vk+i, 

2Ck,  for  rfc+i  <  r  <  r^-i, 

^  /  \  _  /  0’  for  ri  <  r  <  rjv_2, 

j  for  rAr_2  <  r  <  rjv-i, 

where 

Hk{x)  -hk  I  fk{x)dx. 

Jo 

For  the  quadratic  interpolation  function  (2.3.9),  we  get 


Hkix)  =  i2.2Ck+i  -  Ck)x^  +  {Ck  -  Ck+i)x\  k  =  1,2,  ...,JV  -  2;  (2.3.15) 

whereas  for  the  positive- definite  interpolant  (2.3.11),  we  have 

-ff/b(®)  =  Jb  =  l,2,iNr-2.  (2.3.16) 

With  “half  one  one  ...  one  half’  intervals  for  Cjt’s, 

2Ci  =  C2  —  C3  =  ...  =  Cjsf—2  —  SOjv— 1?  (2.3.17) 

we  find  that,  all  Cjk’s  and  hjt’s  are  equal 

Ck  =  2^^  ~  2^^’  ^  ~  l,2,...,iV  —  1.  (2.3.18) 

In  this  special  case,  for  the  quadratic  interpolant  (2.3.9),  we  have  Ak  =  1,  whereas  for  the 
positive  definite  interpolant  (2.3.11),  we  get  m*  =  1.  Thus,  for  both  of  them,  we  obtain 

Hk{x)  =  Ck+ix^  =  ^hkx\  (2.3.19) 


and  recover  the  interpolation  and  blending  functions  reviewed  in  Section  2.2.  In  this  special 
case,  Eqs.  (2.2.6)  and  (2.3.1)  give  the  same  relation  between  C^’s  and  h/t’s.  In  general, 
Eqs.  (2.3.1)  are  different  from  Eqs.  (2.2.6).  In  Fig.  2.3,  we  present  an  example  of  relation 
between  hk's  and  Cfc’s  using  Eqs.  (2.2.6)  and  (2.3.1). 

To  summary  the  result  of  this  section,  we  stress  that  the  scheme  used  here  proceeds 
in  the  following  order: 


(2.3.13c) 

(2.3.14) 


{Ck,  bk)-*{rk,hk)-*{rl^k{r),  Gk{r)). 


(2.3.20) 
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2.4.  Previous  Results  on  C^-Continuity  Blending  Function 

For  blending  functions  Gfc(r)  with  continuity,  the  corresponding  interpolation  func¬ 
tions  are  of  C'  continuity.  It  has  been  shown  in  Ref.  [2]  that,  in  order  to  admit 

the  possibility  of  uniformity  and  to  avoid  any  unspecified  flat  spots,  a  local  interpolation 
function  which  is  not  close  to  either  boundary  point  must  be  non-zero  over  a  minimum  of 
4  intervals.  In  other  words,  the  interpolation  function  V’fc(r)  is  non-zero  in  the  region  of 
^k-2  <  r  <  rk+2  (fc  =  3,4,...,JV’  —  3).  In  Fig.  2.4  the  general  form  of  the  interpolation 
functions  V’t  are  displayed.  As  a  result,  in  the  region  Tk  <  r  <  rjt+i,  only  (up  to)  four 
blending  functions  are  non- zero,  and  the  uniformity  condition  (2.1.12)  reduces  to 


I  ^1  1  ^2  b3  ,  b4 

II  1  J  1 

.s  ,  1 

^1  1 

I  (2.2.6) 

^7 

Cl  Co  C.  Cc 

L‘l  ^  1  M  '  1 

^6  ^7 

1  1  I 

''  1 

>=8 

1  (2.3.1) 

1  "l  1  *>2  1  "3  1  ‘’4  1 

”5  1  he  1 

^1 

^7 

Figure  2.3.  Relations  between  hk's  and  Cfe’s  for  a  case  of  N  =  8. 


'/’k-i{r)+rlfk{r)+ipk+i{r)+i/;k+2{r)  =  1,  r*  <  r  <  rfc+i,  A:  =  2,3,  ...,iV  -  3,  (2.4.1a) 

V'i(r)  +  V'2(r)  +  i>3ir)  =  1,  Ti  <T  <  T2,  (2.4.1b) 

V';v-3(t-)  +  V'N-2(r)  +  V'N-i(’’)  =  1,  rN-2<r<rs-i.  (2.4.1c) 

It  is  easier  to  construct  the  interpolation  functions  starting  from  their  first-order  deriva¬ 
tive,  since  piecewise  continuous  functions  can  be  used.  As  mentioned  in  Section  2.1,  the 
interpolation  functions  must  satisfy  Eq.  (2.1.12).  Consequently,  the  first-order  derivative 
of  all  interpolation  functions  adds  up  to  zero  exactly  in  the  whole  region. 
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(2.4.2) 


N-l 

k=l 

In  Ref.  [2]  V'fc(rfc)  =  0  for  Jb  =  2, 3,  ...,7V  —  2  have  been  used.  Upon  the  evaluation  at  the 
partition  points,  the  uniformity  condition  was  found  to  become 


V’i(n)  +  =  0, 

(2.4.3a) 

V’N-2(^A^-i)  +  V’iV-l(’^7V-l)  =  0, 

(2.4.3b) 

for  end  points  amd 

i>k-ii-rk)  +  i{’k+i(rk)  =  0, 

(2.4.3c) 

for  jb  =  2, 3,  ...,7V -3. 

In  Ref.  [2],  the  same  scheme  as  described  in  Section  2.1  was  used,  i.e.,  the  scheme 
(2.2.12)  was  used  to  first  construct  the  local  interpolation  and  blending  functions  and  then 
find  the  parametric  intervals  Ck  and  locations  bk-  Specializing  to  the  simple  case  of  equal 
intervals  =  fc  (i.e.,  hk  =  1)  and  ak  =  bk  =  Ck  =  dk  =  the  blending  functions  in 
Ref.  [2]  should  reduce  to 


Gkjr) 

Gkir^-i) 


’  n(r  —  k)/24, 
f2(r  -  fc)/25, 

<  [n(r  -  jb)  +  l]/25, 

[n(r  -  ib)  +  f2(r  -  7V)]/11, 
.  [f2(r  -  jb)  +  n(r)  -  37]/ll, 


for  k  =  2,3, ...,  TV  —  2, 
for  k  =  N  —  2, 
for  fc  =  2, 
for  fc  =  TV  —  1, 
for  fc  =  1, 


(2.4.4) 


where 


0, 

for  X  <  —2, 

-2(x  +  2)^ 

for  —2  <*<—§, 

6(x  +  l)’+6(x  + 1)2-1, 

for  — 1  <  X  <  — 1, 

10(x  +  l)2  +6(x  +  l)2  -1, 

for  — 1  <  ® 

-14x2  +  24x  +  12, 

for  —  1  <  X  <  |, 

10(x  -  1)2  -  6(x  -  1)2+25, 

for  1  <  a:  <  1, 

6(x-1)2  -6(x-1)2  +  25, 

for  1  <  X  <  |, 

-2(x  -  2)2  +  24, 

for  1  <  X  <  2, 

24, 

for  X  >  2. 

These  blending  functions  are  illustrated  in  Fig.  2.5.  In  order  to  satisfy  Eq.  (2.1.14),  it  is 
straightforward  to  check  that  the  intervals  must  be 


C3  =  C4  =  ...  =  Cn-3  =  !■> 


(2.4.6a) 
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(2.4.6b) 


Cl  =  Cn-i 


U 
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C2  =  Cn~2  = 


which  are  not  exactly  of  the  “half  one  one  ...  one  half’  spacing  given  in  Section  2.2  for 
blending  function. 


2.5.  New  Results  on  C^-Continuity  Blending  Functions 

Similar  to  the  situation  in  Section  2.3,  questions  here  are:  For  an  arbitrarily  given 
set  of  parametric  locations  bk,  how  to  determine  (or  choose)  the  values  of  hk,  and  how  to 
construct  the  interpolation  functions.  It  is  obvious  that  there  are  many  ways  to  accomplish 
these  things.  But  we  want  to  them  having  certain  properties.  As  in  Section  2.3,  we  want 
to  have  Gk{rk)  =  ^Ck  for  blending  functions  not  close  to  end  points.  We  also  want  to  the 
ratios  of  4  areas  under  the  interpolation  function  to  be  independent  of  the  index  k 

and  to  be  77  :  1  :  1  : 77.  As  will  be  evident  in  the  following  [cf.  Eq.  (2.5.16)],  the  value  of  77 
must  be  within  — 1  <  77  <  0.  [For  blending  functions  (2.4.5),  77  =  —1/13.]  Based  on  these 
considerations,  for  a  given  set  of  Cfc,  we  first  introduce 


Ck 

2(1 +vy 

/:  =  3,4, ...,  N  —  3, 

(2.5.1a) 

^  _  Ck 

k  =  2,N  -2, 

(2.5.1b) 

‘  1  +  277’ 

k  =  1,N  -  1, 

(2.5.1c) 

Co  =  Cu 

Cn  =  Cn-i. 

(2.5.1d) 

Then,  we  choose 

hk  =  Cfc  +  Ck+i  +  +  Ck+2)^  fc  =  1, 2, ...,  N  —2.  (2.5.2) 

The  choice  of  77  should  lead  to  all  hk  positive,  which  requires  that,  for  all  k,  I77I  <  {Ck  + 

Gk+i)/{Ck-i  +  Cfc+2)i  ho-i 

|77|  <  Min{(Cfc  +  Ck+i)/{Ck-i  +  Ck+2)).  (2.5.3) 

In  order  to  write  down  the  area  requirements  in  terms  of  a  unified  form,  we  let 


tpkir)  =  i>k{r),  =  2, 3, ...,  N  -2 


(2.5.4a) 


^i(r)  =  ■0i(r)  +  ^o(’’),  =  rl’N-iir)  +  ^jv(t’).  (2.5.4b) 

After  this,  we  can  rewrite  Eqs.  (2.4.1) 

i>k-i{r) +i>kir) +i>k+i{r) -\-ifk+2{r)  =1,  r*  <  r  <  rt+i,  fc  =  1,2,  ...,iV  -  2.  (2.5.5a) 
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Figure  2.6  shows  the  situation  in  the  interval  <  r  <_rfc+i.  We  also  let  = 

^Nirisf-i)  =  0  so  that  Eq.  (2.1.6)  also  holds  for  xpk,  i.e.,  tAfc(r,)  =  Sjk.  Furthermore,  we 

let 

(n )  =  )  =  0,  (2.5.6) 

as  rl^kirk)  =  0  for  lb  =  2,3,  ...,W  -  2  used  in  Section  2.4.  In  this  way,  we  obtain  from  Eqs. 

(2.4.3) 

i>k-i(rk)  +  ^k+iirk)  =  0,  fc  =  1,2,...,  W  -  1,  (2.5.7) 


Figure  2.6.  Illustration  for  4  non-zero  interpolants  in  the  interval  r*  <  r  <  rjt+i,  k  = 
l,2,...,W-2. 


Corresponding  to  our  choice  (2.5.2),  we  require,  within  the  interval  r*  <  r  <  rk+i,  the  4 
areas  under  the  4  nonzero  interpolants  to  be  determined  by 


rpk-i(T)dr  =  rjCk-i, 


(2.5.8a) 


(2.5.8b) 
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(2.5.8c) 


i>k+iir)dr  =  Ck+i, 


i>k-i-2(r)dr  =  r)Ck+2 


(2.5.8d) 


for  k  =  1,2,  ...,iV  —  2.  It  is  easy  to  see  that,  by  summing  over  both  sides  of  Eqs.  (2.5.8) 
and  using  Eqs.  (2.5.2)  and  (2.4.1),  we  get  an  identity  rt+i  —  r*  =  hk,  which  is  just  Eq. 
(2.1.13).  We  emphasize  that  the  introduction  of  all  tilted  quantities  is  to  eliminate  the 
special  situation  associated  with  the  two  end  points  ri  and  rjv-i  [see  Eqs.  (2.4.3)]  and  the 
two  end  intervals  ri  <  r  <  r2  and  rN-2  <  r  <  r^-i  [see  Eqs.  (2.4.1)].  It  is  easy  to  see 
that  q  =  0  is  a  special  case,  in  which  the  C^-continuity  blending  functions  reduce  to  the 
-continuity  blending  functions  discussed  in  Section  2.3. 

The  piecewise  linear  and  continuous  V'ib(^)  used  in  Ref.  [2]  cannot  be  used  here,  since 
there  would  be  5  parameters  one  can  choose  but  there  are  6  conditions  [three  from  Eqs. 
(2.5.12)  and  the  other  three  will  be  mentioned  below].  In  our  construction,  we  eliminate 
the  intermediate  point  Wk  and  thus  reduce  the  number  of  regions  to  be  considered.  We 
use  piecewise  cubic  polynomials  for  the  first  order  derivative  of  the  interpolation  functions 
^k{r), 


i’ki^)  =  ^kXk  +  Bkx\, 

(2.5.9a) 

^k-ii^)  =  -Zk  +  Xfc  +  B^xl  +  {Zk  -A^-  B^)xl, 

(2.5.9b) 

+  Atxk  +  Btxl  -  {Zk  +  A+  +  Bt)xl, 

(2.5.9c) 

for  k  =  1,2,  ...,N  —  2.  The  fourth  one  V'jfe+2(^)  determined  by  Eq.  (2.4.2)  and  satisfies 
V’lfc+2(*'fc)  =  There  are  6  unknowns  for  a  given  interval:  three  A’s  and  three  B’s.  In  Eqs. 
(2.5.9),  using  Eq.  (2.4.3c),  we  have  defined 

Zk  —  '0fe+l(^fc)  —  '^k—li‘^k)y  k  —  l,2,...,i\r  —  1, 

(2.5.10) 

which 

obtain 

are  positive,  Zk  >  0.  Carrying  out  integration  once  and  using  Eqs. 

(2.1.6),  we 

i>k{r)  =  1  +  hki^Akxl  +  ^Bkxl), 

(2.5.11a) 

i>k-iir)  =  hk[-ZkXk  -h  \a^xI  +  \bZx\  -  i(Afc  +  B^  -  Zk)xl], 

(2.5.11b) 

*+.(r)  =  ht[  Ztx,  +  ^A+xl  +  +  Bt+  ZOxtl, 

(2.5.11c) 
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CHAPTER  III.  CONTROL  POINT  FORM  OF  ALGEBRAIC  GRID  GEN¬ 
ERATION  WITH  CLUSTERING 

The  multisurface  transformation  reported  in  Chapter  II  can  be  used  to  generate  grids 
algebraictdly,  which  is  called  the  control  point  form  (CPF)  of  ^llgebradc  grid  generation 
[3].  In  this  chapter,  we  report  our  progress  made  in  this  area.  This  chapter  is  organized 
as  follows.  In  Section  3.1,  we  identify  the  problem  with  an  implementation  of  the  control 
point  form  of  algebraic  grid  generation.  In  Section  3.2,  we  discuss  how  to  generate  a  linear- 
increment  straight  line  when  all  control  points  Pfc  are  given  on  a  straight  line.  In  Section  3.3 
we  discuss,  when  applying  the  multisurface  transformation  to  a  2D  grid  generation  within 
a  parallelogram,  how  to  make  the  2D  transformation  bilinear.  In  steady  of  a  parallelogram, 
in  Section  3.4  we  discuss  how  to  get  a  bilinear  transformation  within  a  quadrilateral.  In 
Section  3.5  we  discuss  how  to  deal  with  general  case  of  2D  grid  with  clusterings. 

3.1.  Control  Point  Form  of  Algebraic  Grid  Generation 

Grid  generation  in  a  two-dimensional  space  can  be  stated  as  finding  the  relation  P(^ ,  tj) 
between  the  coordinates  P  =  {x,y)  in  the  2D  physical  space  and  the  corresponding  para¬ 
metric  values  (^,7/)  in  a  rectangle  ^min  <  ^  <  ^max,  Vmm  <  V  < 

When  the  multisurface  transformation  is  used  for  algebraic  grid  generation  in  a  2D 
physical  space,  an  array  of  control  points  P*  is  replaced  by  a  net  of  control  points  Qij 
{i  =  1,2,...,/  and  j  =  1,2,...,/).  Let  the  blending  functions  in  Eq.  (2.1.15)  be  0^(0 
(i  =  1,2, ...,/)  for  the  (  direction  and  U  =  15  2, ...,/)  for  the  t]  direction  (in  a 

parametric  space).  Let  the  boundary  specification  be  given  by 

P«i,>l)  =  PUt-u’l)  =  'rHri)-  (3-1.1) 

Starting  from  the  transfinite  interpolation  (also  called  Coons  patches  or  the  Boolean  sum), 
it  htis  been  found  that  the  control  point  form  [3]  of  algebraic  grid  generation  can  produce 
a  2D  grid  in  the  following  equivalent  way:  First,  we  construct  a  tensor  product, 

J  J 

i=l  j=l 

Then,  we  add  edge  adjustments  for  all  four  boundary  curves. 


P((,v)  =  T({,,)  +  a,(0[v‘(.,)  -  T«,,,))  +  a,(OK(,)  - 

+  A(>?)|u‘(0  -  TU,.;.))  +  l3j{v)lnH0  -  (3.1.3) 

The  way  of  assigning  the  coordinates  of  grid  points  in  the  parametric  space  (com¬ 
putational  domain)  is  very  import  for  regenerating  a  given  grid  using  the  CPF  method. 
In  particular,  in  the  presence  of  clustering,  care  must  be  taken.  Since  a  grid  point  car¬ 
ries  index  (i,j),  it  is  uaually  assume  that  the  coordinate  in  the  parmetric  space  is  simply 
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{^,T))  oc  In  other  words,  the  mesh  in  the  parametric  space  is  a  mesh  with  equal 

spacing.  (Sometimes,  the  computational  domain  is  normalized  to  a  unit  square  so  that 
the  spacing  in  the  i  (or  4)  direction  is  different  from  that  in  the  j  (or  t])  direction.  Fig¬ 
ures  3. 1-3.3  show  three  examples  of  regenerating  2D  grids,  including  given  grids,  the  mesh 
distributions  in  the  computational  domain,  and  the  regenerated  grids  using  Eq.  (3.1.3). 
We  see  that  even  in  the  simple  case  of  a  quadrilateral  (Fig.  3.1),  the  given  grid  is  not 
reproduced.  In  the  following  sections,  we  show  how  to  overcome  this  problem. 

3.2.  Generation  of  a  Linear-Increment  Straight  Line 

As  a  preparation  for  later  sections  of  this  chapter,  we  study  in  this  section  how  to 
generate  a  straight  line  P(r)  that  increases  linearly  with  the  parameter  r.  In  other  words, 
we  want  to  know  how  to  choose  the  control  points  Pjt  and  the  blending  functions  Gk{r)  so 
that  the  curve  P(r)  in  Eq.  (2.1.1)  behaves  like  a  straight  line  and  is  also  a  linear  function 
of  the  parameter  r. 


P(r)  =  Pi-h(r-ri)n,  (3.2.1) 

where  n  denotes  the  direction  of  the  straight  hne.  It  is  easy  to  see  that  all  control  points 
must  be  on  the  same  straight  line. 

Pit  =  Pi +(Pfc-Pi)n,  fc  =  2,3,...,7V.  (3.2.2) 

If  the  values  of  pk's  are  not  given,  we  can,  say,  (i)  use  “half  one  one  ...  one  half’  spacing 
rule  in  physical  space  to  determine  pjt’s  and  also  in  parametric  space  to  determine  (Ck,bk) 
and  (ii)  use  “hj  =  /12  =  ...  =  hs-i"  to  determine  blending  functions  Gfc(r)  reported  in 
Sections  2.2  and  2.3.  If,  on  the  other  hand,  we  are  given  the  values  of  p^,  then  the  question 
becomes  how  to  choose  bk-  We  find  that,  if  we  choose 

bk  =  Pki  I*  =  1,2, ...,  AT.  (3.2.3) 

then,  using  Eqs.  (2.1.1),  (2.1.7)  and  (2.1.8),  we  obtain  the  linear  relation  (3.2.1). 

Using  the  linear-increment  straight  line  (3.2.1),  for  any  given  set  of  (curve)  points 

Ci  =  Pi  +  CjH,  i  =  1,2,3,...  (3.2.4) 

(not  control  points),  we  can  reproduce  these  points  Ci  exactly  by  letting  their  parametric 
values  to  be  -|-ri.  If  we  use  a  set  of  control  points  specified  by  Eq.  (3.2.2)  in  the  physical 
space  and  by  Eq.  (3.2.3)  in  the  parametric  space,  we  can  also  regenerate  (curve)  points  Cj 
exactly  using  the  multisurface  transformation  (2.1.1).  The  advantage  of  using  Eq.  (2.1.1) 
is  that,  by  reproducing  initially,  we  can  subsequently  move  one  or  more  control  points 
Pit  and  thtis  modify  the  po.sitions  of  (curve)  points  Ci,  i  =  1,2,3,...  . 

3.3.  Bilinear  Transformation  for  Grid  Generation  in  a  Parallelogram 

In  this  section,  we  consider  the  case  in  which  the  physical  region  where  a  2D  grid  is  to 
be  built  is  a  parallelogram,  i.e.,  the  physical  region  is  a  four-sided  2D  area  whose  opposite 
sides  are  parallel  and  equal.  In  such  a  case,  the  net  of  control  points  can  be  specified  easily 


26 


(a) 


(b) 


(c) 

Figure  3.1.  (a)  A  given  grid  with  clustering  along  two  edges  within  a  quadrilateral,  (b) 
The  uniform  distribution  of  grid  points  in  the  computational  domain  (i.e.,  a  parametric 
space)  used  in  the  CPF.  (c)  The  grid  regenerated  using  the  control  point  form  of 
algebraic  grid  generation  and  the  uniform  grid  distribution  in  the  computational  domain 
shown  in  (b).  The  sparse  and  thick  net  is  the  net  of  control  points,  and  its  distribution 
in  the  parametric  space  is  of  the  “half  one  one  ...  one  half”  type  given  in  Eq.  (2.3.17). 
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(a) 


(b) 


(c ) 

Figure  3.2.  (a)  An  initial  grid  with  clustering  along  one  edge,  (b)  The  uniform  distribution 
of  grid  points  in  the  computational  domain  used  in  the  CPF.  (c)  The  grid  regenerated 
using  the  CPF  method  and  the  uniform  grid  distribution  in  the  computational  domain 
shown  in  (b).  The  sparse  and  thick  net  is  the  control  net,  and  its  distribution  in  the 
parametric  space  is  of  the  “half  one  one  ...  one  hair*  type  given  in  Eq.  (2.3.17). 
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Figure  3.3.  (a)  A  sheet  of  space  shuttle  grid  (near  the  tail  of  the  space  shuttle,  only  half 
grid  is  shown),  which  has  clustering  along  the  body  of  the  space  shuttle. 
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Figure  3.3.  (b)  The  uniform  distribution  of  grid  points  in  the  parametric  space  used  in 
the  CPF. 
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Figure  3.3.  (c)  The  control  net  attached  to  the  2D  grid  given  in  (a).  Its  distribution  in 
the  parametric  space  is  of  the  “half  one  one  ...  one  halF  type  given  in  Eq.  (2.3.17). 
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Figure  3.3.  (d)  The  space  shuttle  grid  regenerated  using  the  CPF  method  and  the  uniform 
grid  distribution  in  the  parametric  space  shown  in  (b). 
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Qij  =  P(6,»?i)  +  Ui  '6)ei  +(Tjj  -m)c2,  (3.3.1) 

where  Cj  and  62  give  the  two  (non-parallel  in  general)  directions  (non-orthogonal  in  gen¬ 
eral)  of  the  parallelogram,  and  is  the  coordinate  of  the  control  point  Qjj  in  the 

parametric  space.  Now,  similar  to  the  linear  relation  given  in  Eq.  (3.2.1),  we  want  the  2D 
transformation  P(^,t7)  to  be  bilinear  in  ^  and  77, 

P(^,^)  =  P(6,m)  +  (^  -  6)ei  +iv-  »7i)e2.  (3.3.2) 

Let  us  see  how  we  can  get  Eq.  (3.3.2).  First,  we  notice  from  Eq.  (2.1.17)  that,  as  a  general 
r\ile  of  surface  generation,  we  have 

/  J 

t=i  j=i 

Then,  from  the  uniformity  condition  (2.1.18),  we  have 

/  J 

Y,  =  (<  Y  =  V  (3.3.4) 

»=i  i=i 

in  the  2D  case.  Substituting  Eq.  (3.3.1)  into  Eq.  (3.1.2)  and  using  Eqs.  (3.3.3)  ,  we  obtain 

I  J 

T(^,77)  =  P(^i,77i) -f  ^ai(0^iei  -  6^1  +  ^/?j(»7)7?j 62  -  7/162.  (3.3.5) 

t=i  i=i 

With  the  choice  and  5'  =  rjj  for  constructing  the  blending  functions  and  using  Eq. 

(3.3.4),  we  simplify  Eq.  (3.3.5)  to 

T(^,77)  =  +  (^  -  6)ei  +  (77  -  7/1)62.  (3.3.6) 

The  bilinear  relation  (3.3.2)  can  be  obtained  either  by  letting  P(^,77)  =  T(^,77)  directly 
[i.e.,  without  edge  boundary  adjustment  terms  in  (3.1.3)]  or  by  parametizing  the  four 
boundary  lines  according  to  the  linear  increment  rule: 


u'  (0  =  P(6 ,  )  +  (^  -  6  )ei ,  U^(0  =  P(6 ,  r? J-I )  +  (^  -  6  )ei , 

v^(7/)  =  P(4i,77i)  +  (7/ -  771)62,  V^(q)  =  P(0-1, ’ll )  +  (»?- 'f?i)e2.  (3.3.7) 

With  the  CPF  constructed  in  this  way,  we  have  a  bilinear  CPF  transformation  within 
a  parallelogram.  Consequently,  for  any  given  grid  within  a  parallelogram,  we  can  regen¬ 
erate  the  given  grid  exactly  provided  that  we  determine  the  coordinate  (^,7/)  of  each  grid 
(control)  point  in  the  parametric  space  using  Eq.  (3.3.2)  [Eq.  (3.3.1)].  Afterwards,  we  can 
move  one  or  more  control  points  Qij  to  manipulate  the  2D  grid. 
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3.4.  Bilinear  Transformation  for  Grid  Generation  in  a  Quadrilateral 

In  this  section,  we  consider  2D  grid  generation  within  a  quadrilateral,  i.e.,  the  physical 
region  is  a  four-sided  2D  area  whose  boundaries  are  4  straight  lines.  In  such  a  case,  we  first 
have  to  normalize  the  2D  parametric  space  to  a  unit  square,  0  =  <  ^  <  =  1  and 

0  =  T/i  <  q  <  T}j-i  =  1.  The  net  of  control  points  should  be  established  now  according  to 


Qij  =  (1  -^0(1  +  ^i(i  - 

+  (1  -  (3.4.1) 

where  is  the  coordinate  of  the  control  point  Qij  within  the  unit  square.  Now, 

instead  of  the  bilinear  relation  in  Eq.  (3.3.2),  we  want  the  2D  transformation  P(^,7/)  to 
be  in  a  different  bilinear  form. 


=  (i-0(i-^)P(6,^i)+^(i-»7)P(^/-i,i?i)+(i-0»fP(6,»/j-i)+^T/P(^/-i,»7J-i); 

(3.4.2) 

i.e.,  a  bilinear  relation  based  on  4  corner  points.  Let  us  see  how  we  can  get  Eq.  (3.4.2). 
Substituting  Eq.  (3.4.1)  into  Eq.  (3.1.2),  we  obtain 


m.v) 


P((uVl) 


«i(0(i  -  (i) 


j 


J=1 


(3.4.3) 


With  the  choice  bi  =  (i  and  6'-  =  rjj  for  constructing  the  blending  functions  and  using  Eqs. 

(3.3.3)  and  (3.3.4),  we  arrive  at 


(3.4.4) 

Similar  to  the  situation  in  Section  3.3,  the  bilinear  relation  (3.4.2)  can  be  obtciined  either 
by  letting  P(^,t7)  =  T(^,t/)  directly  or  by  parametizing  the  four  boundary  lines  as 

u^(0  =  (1  -  <f)P(6,7i)  +  ^P(^/-i,»7i),  u"(0  =  (1  -  0P(6,^7-i) 

v^(t?)  =  (1  -v)^Ui,Vi)  +v'PUi^VJ-i),  =  (1  -  t/)P(^/-i,t7i) -f- T7P(0_i,T?j_i). 

(3.4.5) 
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Having  a  2D  CPF  relation  built  in  this  scheme,  we  get  a  bilinear  CPF  transformation 
within  a  quadrilateral.  Thus,  for  any  given  2D  grid  within  a  quadrilateral,  we  can  regener¬ 
ate  the  given  grid  exactly  provided  that  we  solve  the  coordinate  (^,77)  of  each  grid  (control) 
point  in  the  parametric  space  using  Eq.  (3.4.2)  [Eq.  (3.4.1)].  The  grid  manipulation  can 
be  achieved  by  moving  one  or  more  control  points  Qij  at  a  user’s  discretion.  Note  that 
the  results  in  this  section  can  be  regard  as  a  generalization  of  the  results  in  Section  3.3, 
since  a  parallelogram  is  in  fact  a  special  quadrilateral. 

For  a  2D  grid  with  clustering  either  within  a  parallelogram  or  within  a  quadrilateral, 
it  is  easy  to  see  from  either  Eq.  (3.3.2)  or  Eq.  (3.4.2)  that  the  image  of  the  grid  in  the 
parametric  space  is  also  a  grid  with  clustering.  This  result  shows  that,  for  a  grid 

with  clustering,  its  image  in  the  parametric  space  (^,7?)  should  not  be  with  equal  spacings; 
i.e.,  the  mesh  in  the  parametric  space  should  not  be  like  those  displayed  in  Figs.  3.1(b), 
3.2(b),  and  3.3(b).  In  other  words,  we  should  use  three  spaces  in  the  control  point  form 
of  algebraic  grid  generation.  In  stead  of  the  usual  two  spaces  [physical  space  {x,y)  and 
curvilinear  space  (^,7/)  oc  (t,i)],  we  should  use  the  following  three  spaces:  the  physicad 
region  {x,y)  (arbitrary  shape),  a  parametric  space  (a  rectangle),  and  an  index  space 
(i,jf)  (also  a  rectangle).  In  the  index  space  (t,y),  the  grid  points  are  always  uniformly 
distributed  and  the  grid  lines  are  always  straight  lines.  In  the  parametric  space  (^,77), 
depending  on  the  situation  in  the  physical  region,  the  grid  points  can  be  either  uniformly 
or  non-uniformly  distributed  and  the  grid  lines  may  or  may  not  be  straight  lines.  In  Fig. 

3.4  we  show  how  to  regenerate  exactly  the  grid  given  in  Fig.  3.1(a),  which  is  within  a 
quadrilateral. 

3.5.  General  2D  Physical  Regions 

The  conclusion  reached  at  the  end  of  Section  3.4  is  valid  not  only  when  the  physical 
region  is  a  quadrilateral  (it  contains  the  parallelogram  as  a  special  case)  but  also  for  a 
physical  region  of  any  shape;  Namely,  we  need  three  spaces  for  the  CPF  to  regenerate  a 
given  grid. 

For  a  given  2D  grid  built  within  a  general  2D  region,  we  cannot  have  a  bilinear 
transformation  like  those  in  Eqs.  (3.3.2)  and  (3.4.2).  In  order  to  regenerate  the  given  grid 
using  the  control  point  form  of  algebraic  grid  generation,  we  should  first  reproduce  the 
clustering  patterns  of  the  given  grid  in  the  parametric  space  (^,77)  using,  say,  arc  length 
measurements.  In  general,  such  an  arc-length  method  will  not  reproduce  a  given  grid 
exactly.  It  will,  however,  capture  the  general  clustering  tendency  of  the  given  grid.  Figure 

3.5  shows  how  to  regenerate  the  grid  given  in  Fig.  3.2(a).  In  the  region  of  90-degree  turn, 
the  CPF  grid  does  not  match  the  given  grid  exactly;  Otherwise,  the  CPF  grid  matches 
the  given  grid  exactly.  Figiire  3.6  shows  how  to  regenerate  the  grid  given  in  Fig.  3.3(a)  (a 
sheet  of  grid  for  the  space  shuttle).  Comparing  Fig.  3.6(b)  with  Fig.  3.3(d),  we  see  the 
improvements  near  the  convex  regions  along  the  edge  with  clustering  (i.e.,  77  =  77inin  edge). 
To  reproduce  an  arbitrarily  given  grid  exactly,  we  have  to  subsequently  search  and  adjust 
the  parametric  coordinate  (^,77)  of  each  grid  point  so  that  the  gird  lines  in  the  parametric 
space  are  not  straight  lines  in  general. 
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(b) 


Figure  3.4.  Regenerating  the  grid  of  Fig.  3.1(a)  using  an  improved  implementation  of 
the  CPF  method,  (a)  The  non-uniform  grid  distribution  in  the  parametric  space 
determined  from  Eq.  (3.4.2).  (b)  The  grid  regenerated  using  the  control  point  form  of 
algebraic  grid  generation  and  the  non-uniform  grid  distribution  in  the  parametric  space 
shown  in  (a).  The  sparse  and  thick  net  is  the  control  net,  and  its  distribution  in  the 
parametric  space  is  of  the  “half  one  one  ...  one  halF  type  given  in  Eq.  (2.3.17). 
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(a) 


(b) 


Figure  3.5.  Regenerating  the  grid  of  Fig.  3.2(a)  using  an  improved  implementation  of 
the  CPF  method,  (a)  The  non-uniform  grid  distribution  in  the  parametric  space 
determined  according  to  the  grid  spacing  in  the  physical  space,  (b)  The  CPF  grid 
generated  using  the  non-uniform  grid  distribution  in  the  parametric  space  shown  in 
(a).  The  sparse  and  thick  net  is  the  control  net,  and  its  distribution  in  the  parametric 
space  is  of  the  “half  one  one  ...  one  half’  type  given  in  Eq.  (2.3.17). 
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Figure  3.6.  Regenerating  the  grid  of  Fig.  3.3(a)  using  an  improved  implementation  of 
the  CPF  method,  (a)  The  non-uniform  grid  distribution  in  the  parametric  space 
obtained  through  an  arc-length  measurement  of  the  space  shuttle  grid  in  the  physical 
space. 
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Figure  3.0.  (b)  The  space  shuttle  grid  regenerated  using  the  CPF  method  and  the  non- 
uniform  grid  distribution  in  the  parametric  space  shown  in  (a).  The  sparse  and  thick 
net  is  the  control  net,  and  its  distribution  in  the  parametric  space  is  of  the  “half  one 
one  ...  one  half’  type  given  in  Eq.  (2.3.17). 


CHAPTER  IV  CURVATURE  CORRECTION 


4.1  Basic  Concepts  of  Elliptic  Grid  Generation 

The  underlying  idea  behind  all  grid  generation  methods  is  to  express  the  physical 
coordinates  of  a  point  in  space,  in  terms  of  a  set  of  variables  in  the  computational  domain.  For 
each  set  of  physical  coordinates  (x,y,z),  there  is  a  corresponding  set  of  points  in  the 

transformed  or  computational  domain.  The  mathematical  expression  of  the  above  statement  is: 

x  =  x(^,Ti,0  <=>  4  =  ^(x,y,z) 

y  =  y  <=>  Ti  =  Tj  (x,y,z) 

z  =  z(4,-n,0  «  C  =  C(x.y,z) 

Among  the  many  techniques  available  for  numerical  grid  generation,  elliptic  techniques  are 
perhaps  the  most  widely  used.  We  restrict  the  discussion  here  to  two  dimensions.  We  start  by 
writing  a  set  of  elliptic  equations  for  the  transformed  variables  (^."n).  The  simplest  elliptic  equation 
is  Laplace's  equation; 

4xx  +  ^yy  =  0  (4.1.1) 

Tlxx+Tlyy=0  (4.1.2) 

These  equations  state  that  given  values  of  the  transformed  variables  ^  and  tl  or  their  derivatives  on 
the  boundary  of  an  arbitrary  physical  domain,  we  can  compute  values  of  ^  and  T|  inside  the 
domain.  However,  our  goal  is  to  express  x  and  y  inside  the  physical  domain  in  terms  of  ^  and  T|. 
Let  us  assume  the  computational  or  transformed  domain  is  rectangular  and  ^  and  "n  are  distributed 
uniformly  over  this  domain.  Let  the  increments  j  and  At|  =  •ng+i  -  T|ij_i  be 

equal  to  2  for  convenience.  The  computational  or  transformed  domain  is  thus  a  rectangular  region 
with  a  rectangular  mesh  covering  it.  Our  goal  is  to  find  values  for  x  and  y  at  each  point  on  the 
computational  domain,  given  values  of  x  and  y  or  their  derivatives  on  the  boundaries. 

Exchanging  the  role  of  the  dependent  and  independent  variables  in  Equations  (4.1.1)  and 
(4. 1 .2)  transforms  this  set  of  equations  to: 

ax^^  -  2Px^ti  +  =  (^  (4.1.3) 
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(4.1.4) 


where  the  terms: 


a  =  +  yri^ 

(4.1.5) 

P  =  x^Xt^  +  y^y^, 

(4.1.6) 

Y=x^  2^y^2 

(4.1.7) 

are  responsible  for  the  coupling  between  the  two  equations. 


It  should  be  immediately  apparent  that  the  transformed  equations  do  indeed  represent  equations  for 
X  and  y  in  terms  of  4  and  r\. 

Equations  (4.1.3)  and  (4.1.4)  are  discretized  and  solved  simultaneously  for  x  and  y.  The 
discrete  form  of  the  equations  is  simple,  because  the  spacing  in  the  computational  domain  is 
uniform  and  the  increments  are  equal  to  1.  For  example,  the  discretized  form  of  (4.1.3),  based  on 
central  differences  is: 


a. .  = 
‘0 

-X.  .  ,)f 
llj-l 

Py- 

Xi.y)(X. 

,i.i  -  ■  ^i. 

j-1^ 

■  X.  ,  .)f 

a. .  (x.^,  .  -  2x. 

i,j '  i+i,j  1 

i.j'*' Vi.j 

-  ^ 

■Vij+r 

X 

T 

+ 

X 

T 

II 

0 

(4.1.8a) 

Similarly  for  (4.1.4): 

“i,j  -  2yi,j  *  yi-ij )  *  r,,j  (yy*,  -  zyu  +  yy-,) 

-  2P1.J  yi*!,,-!  *  yi-i.M> = ® 


(4.1.8a) 


4.2  Properties  of  the  Elliptic  System 

Elliptic  systems  exhibit  the  desirable  tendency  to  smooth  out  the  grid  insioe  a  domain,  even 
if  the  boundary  is  not  smooth.  They  also  guarantee  that  the  grid  lines  do  not  cross.  These 
properties  lead  to  their  popularity. 
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The  behavior  of  elliptic  systems  near  curved  boundaries  will  be  the  focus  of  our 
discussion.  Elliptic  grids  behave  differendy  near  concave  and  convex  boundaries.  In  general,  they 
behave  well  if  the  boundary  is  convex  when  viewed  from  the  interior  of  the  domain  and  behave 
poorly  if  the  boundary  is  concave.  This  is  because  the  grid  lines  tend  to  cluster  near  convex 
regions  at  the  expense  of  concave  regions.  This  behavior  manifests  itself  in  two  ways:  If  the 
points  on  the  boundaries  are  fixed,  the  normal  spacing  next  to  concave  boundaries  can  becomes 
orders  of  magnitude  larger  than  the  spacing  next  to  convex  boundaries.  If  the  points  are  allowed  to 
float  along  the  boundary,  in  addition  to  the  previous  phenomena,  the  points  tend  to  migrate  away 
from  concave  comers  along  the  boundary,  resulting  in  poor  resolution  of  the  boundary  curve. 

4.3  Curvature  Control 

The  deterioration  of  grid  quality  resulting  from  the  behavior  of  elliptic  systems  at  curved 
boundaries  is  severe  enough  to  warrant  corrective  measures.  The  most  intuitive  approach  is  to  deal 
with  the  elliptic  system  on  a  discrete  level  and  try  to  find  the  terms  responsible  for  this  behavior. 
Once  the  guilty  terms  have  been  identified,  they  can  be  eliminated  or  damped,  keeping  in  mind  that 
eliminating  the  ill  effects  of  these  terms  near  concave  comers,  also  diminishes  their  favorable  effect 
along  convex  comers.  Our  primary  goal  will  be  to  moderate  the  effect  of  curvature  on  the  grid,  not 
to  eliminate  it. 

Obtaining  a  numerical  solution  to  the  pair  of  elliptic  equations  (4.1.8a)  and  (4.1.8b)  begins 
by  starting  with  an  initial  grid  whose  point  distribution  does  not  satisfy  those  equations.  The 
points  of  the  initial  grid  must  move  to  a  new  position,  until  they  satisfy  equations  (4.1.8a)  and 
(4.1.8b)  and  their  boundary  conditions.  We  will  try  to  understand  the  exact  nature  of  the  forces 
that  cause  the  points  which  do  not  satisfy  the  elliptic  equations,  to  move  until  the  grid  does  satisfy 
the  elliptic  equations.  We  will  confine  our  discussion  to  grids  with  orthogonal  boundaries,  and 
restrict  the  regions  where  we  wish  to  exercise  curvature  control  to  the  vicinity  of  the  boundaries. 
For  regions  of  the  grid  which  are  orthogonal; 

(J  =  x^Xt^  +  y^yr,  =  0  (4.3.1) 

Notice  that  although  an  orthogonal  grid  satisfies  equation  (4.3.1),  simply  dropping  P  in  equations 
(4.1.8a)  and  (4.1.8b)  does  not  force  equation  (4.3.1)  to  be  satisfied  and  therefore  will  not 
guarantee  an  orthogonal  grid.  Assuming  that  we  have,  by  some  means,  forced  the  grid  to  be 
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orthogonal  in  a  region,  we  can  then  set  P  =  0  in  that  region.  This  simplifies  equations  (4.1.8a)  and 
(4.1. 8b)  considerably : 


“ij  *  ^i-ij )  +  (yg+i  -  *  ^g-i)  “  ® 

Solving  equation  (4.3.2)  for  x. .  and  dropping  the  i,j  subscript  from  a. .  and  y. 

IJ  ItJ  IJ 


A  similar  expression  can  be  written  for  yj 


For  regions  of  the  grid  which  are  orthogonal,  the  point  at  (i  j)  depends  only  on  its  four  neighbors 
at  (i+1  J),  (i-1  J),  (ij+l)  and  (i,j-l).  As  was  mentioned  earlier,  this  discrete  form  of  the  elliptic 
equation  forms  the  basis  of  our  attempt  to  isolate  the  terms  responsible  for  the  undesirable  behavior 
of  the  elliptic  system  at  concave  boundaries. 


Typically,  for  all  the  points  in  the  field  which  are  in  orthogonal  regions  of  the  grid,  the 
point  at  (ij)  is  forced  to  move  to  a  new  location,  as  to  satisfy  equations  (4.3.3a)  and  (4.3.3b).  The 
tendency  of  the  points  to  move  away  from  concave  comers,  is  part  of  this  motion.  In  order  to 
moderate  this  tendency,  we  must  understand  the  exact  nature  of  this  movement  and  try  to  isolate  its 
different  components.  Denote  the  position  vector  at  the  point  at  (i  j)  by  Pj  j: 


It  is  evident  that  in  order  for  this  point  to  be  in  an  equilibrium  positions  relative  to  its  neighbors  in 
an  orthogonal  region  of  the  grid,  x. .  and  y^ .  must  be  given  respectively  by  equations  (4.3.3a)  and 

(4.3.3b).  The  relative  motion  vector  Rj^  is  defined  as  the  difference  in  the  position  of  a  point 
before  and  after  it  satisfies  equations  (4.3.3a)  and  (4.3.3b).  In  symbolic  form: 


R. .  =  (P. .) 

'.J  '.J  n 


-  (P.  .) 

»  ».J  old 
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For  a  fully  converged  grid  R.  ^  tends  to  zero  every  where  in  the  domain.  To  distinguish  between 
orthogonal  and  non  orthogonal  regions,  we  use  the  subscript  1  and  2.  satisfies  the  elliptic 

equations  everywhere,  whereas  R2-  •  is  restricted  to  orthogonal  regions.  We  have  thus  identified 

the  vector  which  describes  the  direction  and  magnitude  of  the  motion  of  a  point.  We  must  now 
find  its  different  components. 

Some  manipulation  of  equations  (4.3.3a)  and  (4.3.3b)  shows  that  R2.  can  be  expressed 

as  the  sum  of  two  vectors  which  depend  only  on  the  difference  in  the  current  position  of  the  center 
point  and  its  four  neighboring  points.  If  we  let  the  vector  Vj  contain  the  contributions  from  i  +  1 
and  i  - 1  and  the  vector  contain  the  contribution  from  j  +  1  and  j  -  1,  the  following  holds  true: 

R2,-V^-.V,  (4.3.4) 

Figure  4.1  is  a  graphical  representation  of  the  elliptic  equations,  showing  the  contribution  of  the 
neighboring  points  to  the  final  location  of  the  center  point 


Fiturc  4.1 


To  derive  equation  (4.3.4)  formally,  it  is  convenient  to  define  two  new  scalar  quantities  a 
and  b,  such  that: 


a 

a  = - 

a  +  Y 


b  = 


y 


a  +  Y 


44 


using  these  definitions  in  (4.3.3a)  and  (4.3.3b),  and  subtracting  (x. .)  and  (y. .)  from  both 

‘J  old  *0  old 

sides  of  these  equations  respectively,  the  Cartesian  components  of  the  vectors  Vj  and  can  be 
expressed  as: 


Vj  =|a(Ax|-Axj  )i+^a(Ayj-Ay.)j  (4.3.5a) 

=  I  b  (AxJ  -  Ax‘  )  i  + 1  b  (Ay]  -  AyT)  j  (4.3.5b) 


Where  the  following  shorthand  notation  is  used: 


Ax.  =  x.^,  .  -  X. . 

1  1+1. j  i.j 


Ax.  =  X.  ,  .  -  X. . 
1  1-1. J  i.j 


Ax  .  =  X. .  ,  -  X. . 
J  i.J+1  1.J 


Ax.  =  X. .  ,  -  X. . 
J  i.J-1  M 


‘^>'1  =  VKi-yy 

dy  -  =  yy^i  -  yy 
Ay~  =  y. .  -  y. .  , 


Notice  that  in  the  above  definitions,  the  subscript  ’old’  has  been  dropped  from  the  points  at  i,j  to 
make  the  equations  easier  to  read.  Also  keep  in  mind  that  in  Figure  4.1  the  effect  of  the  quantities  a 
and  b  on  Vj  and  V2  have  not  been  taken  into  consideration. 

Adding  equations  (4.3.5a)  and  (4.3.5b)  gives  the  Cartesian  components  of  R2ij' 

R2i.j  =  Rx'  +  RyJ 


where 


Rx  =  |a  (Ax]  -  Ax.  )  + 1 b  (Ax]  -  Ax^  ) 


Ry  =  ^a  (Ay]  -  Ayj  )  +  ^  b  (Ay]  -  Ay^  ) 


(4.3.6a) 

(4.3.6b) 


Although  equations  (4.3.6a)  and  (4.3.6b)  give  the  components  of  point  movement  in  Cartesian 
coordinates,  in  this  form  it  is  not  easy  to  determine  which  of  these  components  pulls  the  point  in  a 
particular  direction  of  the  physical  domain.  As  a  particular  example  consider  a  solid  boundary 
located  along  the  ^  =  constant  direction.  Since  the  ^  direction  runs  along  a  constant  value  of  the  j 
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index,  it  becomes  obvious  which  components  of  the  R2j  j  vector  pull  the  points  away  from  the 
surface  if  the  components  of  the  vector  are  in  the  %  and  "n  or  i  and  j  directions.  (The  i  and  j  indices 

A  A  -» 

do  not  point  in  the  i  and  j  directions.)  An  alternative  expression  for  R2.  •  as  a  sum  of  two  vectors 

•0 

in  the  ^  and  r\  directions  is  therefore  highly  desirable.  Figure  4.2  suggests  clearly  that  the  two 

A  A 

vectors  we  are  seeking  are  a  and  y.  which  do  indeed  point  along  the  T|  and  q  directions 

*’■*  A  A 

respectively.  A  formal  definition  of  the  two  unit  vectors  a. .  and  Y  .  is: 

UJ  IJ 

(4.3.7a) 

(4.3.7b) 


^  >0 

A  ]  A  A 


* 


FIcurc  4.2 


The  quantities  a  ^  and  Yj  ^  are  formally  defined  as: 


a. .  =  (P„  •  P„).  . 

I.J  T1  TT'i.J 


We  can  thus  write: 

R2j  .  =  R,  i  +  Ry  j  =  (Viy  +  V2y)yj .+  (Via  +  V2a)aj  j  (4.3.9) 
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The  quantities  >^la  ^2a  components  of  the  vectors  Vj  and  in  the  and  a 

A  A  A 

and  Y  directions.  For  arbitrary  unit  vectors  a  and  y  which  make  an  angle  9  between  them,  we 
can  write  by  consulting  Figure  4.2: 

ViY=Vj  ■  Y- Videos 0 

Vi„  =  Vj  •  a  -  Vi^yCos  0 

similarly: 

V2y=V2- Y- V2aCOS0 

—>  A 

Via  =  •  a  -  ViyCos  0 

The  above  set  of  equations  can  be  solved  for  the  four  unknowns  Vj^ ,  Viy  .Vja  and  Via- 
However,  we  mentioned  that  our  grid  is  orthogonal  in  the  regions  in  which  we  are  interested.  This 
means  the  angle  0  is  90°.  As  a  result,  the  above  set  of  equations  reduce  to  taking  the  dot  product 
of  V,  and  V,  with  a. .  and  y. ..  Using  equations  (4.3.5a)  and  (4.3.5b)  along  with  (4.3.7a)  and 

(4.3.7b),  we  get: 


In  Figure  4.2,  if  curvature  at  the  point  (i  j)  where  to  have  no  pulling  effect  on  the  grid  point 
at  (i  j),  the  pair  of  points  at  (i+lj)  and  (i-lj)  should  pull  the  central  point  along  the  a. .  direction 

A 

only  and  similarly  the  pair  of  points  at  (iJ+1)  and  (iJ-1)  should  pull  it  along  the  Yj  j  direction  only. 
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In  other  words  by  setting  V  and  ^2-^  equal  to  zero  the  pulling  effect  of  curvature  will  be 
completely  eliminated. 

Remember  however  that  we  are  going  to  be  prudent  and  try  to  retain  as  much  of  the 
curvature  effect  as  possible.  Also,  in  certain  regions  we  might  want  to  diminish  the  effect  of 
curvature  selectively.  In  terms  of  our  equations,  the  above  observation  amounts  to  rewriting 
equation  (4.3.9)  in  the  following  form: 

R2j  .  =  R,  S  +  Rj,  J  =  (V,^  +  (OV,a  +  V2J0,  j  (4.3.14) 

In  this  form,  adjusting  the  values  of  P  and  Q  allows  us  to  adjust  the  pulling  effect  of  curvature  in 
each  direction  independently.  But  we  are  not  finished!  We  have  so  far  worked  under  the 
assumption  that  the  grid  is  orthogonal  in  the  regions  we  are  interested  in.  What  happens  in  other 
regions?  Obviously  as  the  grid  deviates  from  orthogonality,  the  accuracy  of  equations  (4.3.3a)  and 
(4.3.3b),  which  formed  the  basis  of  our  derivation  and  led  to  (4.3.14),  diminishes.  We  need  a 
way  to  compensate  for  this  and  provide  a  smooth  transition  between  the  regions  which  are 
guaranteed  to  be  orthogonal  and  those  which  are  not 

Let  the  relative  motion  vector  R. .  be  composed  of  two  components,  one  which  satisfies 
the  full  elliptic  equations;  that  is  equations  (4.1.8a)  and  (4.1.8b)  and  is  donated  Rjj ^  and  another 
which  satisfies  equation  (4.3.10)  and  is  donated  R2. .,  consistent  with  our  previous  definitions  of 
these  two  vectors.  Let  a  scalar  <I)  be  defined  such  that  it  varies  smoothly  between  o  and  1.  We  can 
write: 

R..  =  <DR,..+  (1 -<D)R2..  (4.3.11) 

1.J  MJ  ' 

Equation  (4.3.1 1)  states  that  unless  <I>  is  equal  to  0,  the  solution  will  not  be  without  the  influence 
of  the  full  elliptic  equations,  and  that  if  C>  is  equal  to  1,  the  solution  depends  entirely  on  Ri  j  j- 

4.4  The  GridPro™/pc2000  Computer  Program 

The  above  technique  of  curvature  control  has  been  successfully  implemented  in  the 
GridPro™/pc2000  elliptic  grid  generation  computer  program.  This  program  forces  the  grid  to  be 
orthogonal  at  the  boundaries  of  the  domain.  In  the  interior  regions  of  the  grid,  the  orthogonality 
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condition  is  not  enforced,  allowing  deviation  from  orthogonality.  In  most  practical  cases  the  grid 
stays  very  close  to  orthogonal  even  in  the  interior  of  the  domain. 

Equation  (4.3.11)  is  implemented  in  GridPro™/pc2000  by  letting  <I>  vary  exponentially 
between  a  minimum  value  at  the  walls  and  1  in  the  interior  regions.  The  minimum  value  depends 
on  the  topology.  For  O  grids  the  minimum  value  is  0.25.  The  scalars  P  and  Q  are  also  adjusted 
according  to  the  topology.  Again  we  consider  the  O  grid  as  an  example.  Here  the  direction  along 
the  body  is  the  ^  and  the  direction  normal  to  the  body  the  T|  direction.  When  an  elliptic  grid  with  an 
O  topology  is  constructed  about  a  body  such  as  an  ellipse,  the  region  near  the  tip  shows  an  intense 
clustering  of  points,  as  seen  in  Figure  4.3a.  This  clustering  is  desirable  in  the  ^  direction,  but  the 
grid  can  be  improved  by  reducing  the  clustering  in  the  T|  direction.  When  this  is  done,  the  result  is 
the  grid  in  Figure  4.3b.  Notice  the  normal  spacing  is  now  much  more  uniform  along  the  ellipse. 
The  grid  in  Figure  4.3a  shows  a  significant  change  in  normal  spacing  between  the  top  part  and  tip 
section  of  the  ellipse.  The  grid  in  Figure  4.3b  was  produced  using  GridPro™/pc2000  by  setting  Q 
equal  to  0.  Figure  3a  was  produced  with  all  the  settings  the  same  except  with  Q  equal  to  1.  For 
both  cases  P  was  equal  to  1. 


Figure  4.3a 


Figure  4.3b 
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CHAPTER  V.  SUMMARY 

During  our  Phase  II  work,  (i)  we  have  made  the  control  point  form  (CPF)  of  algebraic 
grid  generation  more  applicable  in  the  real  world  of  grid  generation,  and  (ii)  we  have 
developed  an  effective  method  to  control  the  grid  point  distribution  of  elliptic  grids  along 
curved  boundaries. 

(i)  Progress  made  in  the  CPF  of  algebraic  grid  generation  was  reported  in  Chapters 
II  and  III.  In  Chapter  II,  we  reviewed  the  old  blending  functions  and  described  how  to 
construct  new  and  more  flexible  blending  functions  as  well  as  the  results,  both  for  C^- 
and  C^-continuity  blending  functions.  The  new  blending  functions  can  allow  a  computer 
software  user  to  choose  the  locations  of  control  points  arbitrarily.  In  Chapter  III,  we 
discussed  how  to  make  the  2D  transformation  relation  given  by  the  CPF  to  be  a  bilinear 
transformation  within  a  quadrilateral  (which  includes  a  parallelogram  as  a  special  case). 
This  helps  us  to  recapture  the  general  clustering  feature  of  a  given  grid  of  arbitrary  shape 
when  we  use  the  CPF  to  regenerate  the  given  grid.  These  two  improvements  have  been 
included  in  one  of  our  grid  generation  programs  called  GridPro/sb,  which  can  be  used  to 
improve  the  quality  of  a  single-block  3D  volumetric  grid.  In  our  software  GridPro/sbSOlO 
and  GridPro/sb3015,  a  user  can  insert  or  remove  a  control  surface  (or  a  control  line  in  a 
boundary  surface)  one  at  a  time.  In  software  GridPro/sb3000-GridPro/sb3015,  a  user 
has  the  choice  of  preserving  the  clustering  of  an  initial  grid. 

(ii)  Progress  made  in  elliptic  grid  generation  was  reported  in  Chapters  IV.  Our  method 
of  controling  the  grid  point  distribution  near  curved  boundaries  has  been  implemented  in 
our  grid  generation  code  GridPro/pc2000.  This  code  produces  more  uniform  grid  point 
spacing  along  curved  boundaries  in  both  the  tangential  and  normal  directions. 
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